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Donald  P.  Gaver 
Patricia  A.  Jacobs 

Naval  Postgraduate  School 
'  Monterey,  California  93943 
U.S.A. 


This  paper  addresses  two  problems  of  Interest  In  service 
system  analysis:  (a)  that  of  making  statistical,  data- 
driven  estimates  of  the  long-run  probability  of  a  long 
delay,  and  (b)  the  assessment  of  rate  of  approach  to  a 
long-run  system  performance  measure  such  as  expected 
delay,  the  rate  being  characterized  by  a  simple  exponen-* 
tial,  at  least  Initially.  Both  ^e  Illustrated  by 
reference  to  M/G/1  and  related  systems. 
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The  application  of  probability  theory  to  a  wide  variety  of*  congestion  problems 
arising  In  communication  systems  has  been  well  catalogued  In  many  papers;  the 
treatises  by  Borovkov  (1984),  Cohen  (1969),  Cooper  (1972),  Cox  and  Smith  (1961), 
Feller  (1966),  Kleinrock  (1975),  Gross  and  Harris  (1974),  and  many  others  attest 
to  the  attractiveness  of  the  area  for  probability  modelling.  In  this  body  of  work, 
several  features  are  noticeable.  First,  most  of  the  elegant  solutions  obtained 
are  In  somewhat  Implicit  form,  being  presented  as  functional  equations,  or,  fre¬ 
quently,  as  Integral  (Laplace)  transforms,  generating  functions,  and  sometimes  as 
combinations  of  the  above.  Moments  (particularly  the  first)  of  delays  or  number 
In  queue  under  long-run  conditions  when  steady  state  prevails  are  the  most  explicit 
and  assessible  performance  summaries.  Second,  the  results  obtained  are  presented 
In  terms  of  component  distribution  functions  and  stochastic  processes  (renewal, 
Poisson,  etc.)  that  are  taken  as  known;  only  rarely  are  Issues  addressed  that 
arise  when  actual  data  Is  to  be  used  as  a  basis  for  Inference  from  the  models; 
however,  see  Cox  (1965).  Third,  only  fragmentary  comprehensible  Information  con¬ 
cerning  transient  behavior,  not  to  mention  time-dependence  of  parameters  Is  avail¬ 
able;  recently  however  work  by  Roth  and  Odoni  (1983),  Roth  (1981),  Abate  and 
Whitt  (1985),  Lee  (1985)  has  elucidated  the  simple  exponential -like  approach  to 
the  steady  state  displayed  by  a  great  many  stochastic  systems  with  stationary, 
time-homogeneous  arrivals  and  service  processes. 


Instead  of  directly  dealing  with  the  above  classes  of  problems  much  attention 
has  been  concentrated  on  modelling  large  networks  of  servers,  particularly  steady- 
state  Markovian,  cf.  Kelly  (1979)  for  an  elegant  treatment;  see  package  programs 
created  by  AT4T  (the  D.  Mitra  PANACEA)  and  by  IBM  (work  by  Lavenberg,  ^  al). 

Other  numerical  work  by  Neuts  (1981)  and  associates  has  pioneered  the  exploita¬ 
tion  of  special  structure  of  certain  stochastic  systems.  And  a  considerable  ef¬ 
fort  to  create  approximations,  ranging  from  heavy-traffic  and  diffusion  to,  lately, 
WKB  approaches,  see  Knessl,  Matkowsky  ^  al_  (1986)  has  been  evident;  see  in  par¬ 
ticular  Newell  (1982)  for  pioneering  work. 

This  paper  deals  specifically  with  approaches  to  the  two  classes  of  somewhat 
neglected  problems  referred  to  earlier:  those  of  data-driven  inference,  and  of 
assessing  transient  behavior.  The  approaches  proposed  are  Illustrated  concretely 


In  terms  of  the  simple  M/G/1  system,  but  apply  more  widely. 

2.  Inference  Concerning  Long  Steady-State  Delays  in  M/6/1 -Like  Systems 

Consider  a  single  service  system  approached  by  stationary  Poisson  (A)  traffic 
with  A  known.  Service  times,  or  message  lengths,  X,  are  Independently  distributed 
as  Fx(x);  assume  AE[X]  =  p  <  1.  This  Is  the  M/G/1  system.  Let  observations  of 
the  service  times  be  all  that  is  known  about  F:  denote  these  by  X] ,X2<* • • tX^. 

The  objective  Is  to  supply  estimates  of  long-run  characteristics  of  the  system: 

In  particular,  estimates,  and  error  assessments  thereof,  of  the  probability  of  a 
long  delay  experienced  by  an  arriving  message. 


2.1  Virtual  Waiting  Time  Transform;  The  Long  Tall. 


It  Is  well  known  that  If  W(t)  Is  the  virtual  waiting  time  In  the  M/G/1  and 
p  <  1,  then  the  Laplace  transform 


ECe'*''] 


llm  • 

t-»<» 


1"P 


se!xT^J 


1-0 

l-pA(s) 


(2.1) 


where  A(s)  is  the  Laplace-Stleitjes  transform  of  a  distribution.  It  can  also  be 
seen  that  In  various  Interrupted  and  priority  service  systems,  including  those  in 
which  the  server  "takes  vacations,"  that  the  above  can  be  changed  to 


£[6"^“]  =  -CCs)  (2.2) 

l-pB(s) 

where  possibly  B(s)  is  the  transform  of  a  completion  time  (effective  service  time, 
i.e.^X  including  durations  of  a  random  number  of  interruptions  occurring  therein), 
and  C(s)  is  the  transform  of  an  honest  distribution  that  accounts  for  effects  of 
busy  period  starts;  see  Gaver  (1968). 


If  A(s)  (or  B(s))  exists  for  s  >  Sq,  Sg  <  0  then  there  will  be  a  smallest 
real  zero,  -s  =  k  >  0,  of  the  denominator  of  (2.1 ),( (2.2) ) ,  which  can  be  used  to 
show  that 


P{W  >w}  ~  D(K)e''^,  w-*-®  (2.3) 

One  way  of  establishing  the  above  exponential  tail  property  is  to  introduce 


'I'  (s) 


1-E[e'^*^] 


=  /  P{W  >w}e"^'^dw 
0 


into  (2.1)  which  leads  to 


'l'(s)  =  p[l~Ms)_]  +  pA(s)ij;(s) 

equivalent  to  the  terminating  renewal  equation,  see  Feller  (1966),  p.  362, 

_  _  _ 

Fy(w)  =  P{W>w}  =  H(w)  +  p  /  F(w-x)A(dx) 

“  0 

Introduce  ^(w)  »  Fy(w)e®'^,  0  real  and  pos'itive,  into  (2.5) 
fJJ(w)  =  H^(w)  +  /  F^(w-x)pe®^A(dx)  , 
and  choose  6  »  <  so  that 


(2.4) 


(2.5) 


(2.6) 


/  pe®^A(dx)  =  /  A*^(dx)  =  1  , 
0  0 


.#/ 


(2.7) 


yielding  a  standard  renewal  equation  for  to  which  the  key  renewal  theorem 
applies,  yielding  as  w  -►  * 


/  r(w)dw 


11m  F?[(w)  =  11m  Fy(w)e'^’^  »  D{K)  »- 
w«>  “  w-^  “  f  J 


(2.8) 


/xA"(dx) 

0 

equivalent  to  (2.3).  A  similar  expression  Is  available  for  (2.2). 

2.3  Statistical  Estimation  of  Probability  of  Long  Walt. 

If  data  are  available  on  service  times  (message  lengths)  then  an  estimate  of 
the  requisite  transform  Is 


-sx 


1 


and  K  Is  the  unique  positive  solution  of  this  sample  equivalent  of  (2.7): 

,  ,  n  0x. 

«  '-!)  =1  • 


(2.9) 


(2.10) 


1=>1 


A  three- term  Taylors'  series  expansion  around  9  »  0  gives  an  Initial  estimate 

'  .  (2.11) 

where  x^  »  (x|f  +  ...  +  x[5)/n,  the  sample  moment}  (2,11)  Is  recognized  as  being 

the  non-parametric  sample  version  of  the  exponential  distribution  parameter  ob¬ 
tained  by  normalizing  to  W/ECW]  and  letting  p  i  1  (the  heavy  traffic  approximation). 
A  search  or  Newton-Raphson  Iteration  starting  from  (2.11)  quickly  yields  <  as 
accurately  as  Is  necessary. 

A 

Some  theoretical  asymptotic  results  concerning  modes  of  behavior  of  k  will  now 
be  sketched;  more  detail  will  be  presented  elsewhere.  First,  It  Is  essential  that 

ECe*^^]  <  00  In  order  for  (2.7)  to  hold,  and  so  by  continuity  f^"’^(<)  =  E[x'’’e'^^]  <oo 
for  any  Integer  m.  Now  express  (2,10)  as 

11^  9x. 

and  expand  In  Taylor''s  series  around  k:  since  f(K:)  »  0,  for  some  (|) 


(2.12) 


0  »  f(K)  +  (K-K)f'(K)  +  -  K)^f "(({)<)  , 


so 


K  -  K  » 


f'(<)  +  ^<-K)f"((j)K) 


(2.13) 

(2.14) 


Note  that  If  E[e^''^]  <  oo  then  VarCe*^^]  <  oo  and 


Var[f(K)]  =  ^  VarCe''^]  , 


(2.15) 


ECf(<)]  *  0  . 

^  Ak  ^ 

Additionally,  f'(K)  -►  E[f’ (k)3  by  the  weak  law  of  large  numbers  and  <-<  -*■  0  In 
probability,  so  by  the  central  limit  theorem, 

-  N(O.l)  (2.16) 

*  K 

2kX 

the  standard  Normal /Gaussian  distribution.  If,  however,  E[e  ]  =  ~,  then  If  any 
kind  of  limiting  distribution  is  to  exist,  the  components 


of  the  sum  f(K)  must  be  in  the  domain  of  attraction  of  a  stable  law,  which  will  be 
assumed.  In  any  case  -►  £[?'(<)]  as  before, so  a  proper  normalization  by  a 

power  of  n  shows  that  (k-k)  has  a  stable  law  as  limiting  distribution.  Fortunate¬ 
ly,  the  relatively  well-behaved  Normal/Gauss  limit  prevails  when  traffic  intensity 
is  relatively  high  (p  >  1/2  in  case  X  ~  Exp(y),  for  example).  ^ 

2.4  Assessing  Variability  of  the  Estimates. 

Having  obtained  an  estimate  of  the  parameter  <,  and  of  the  probability  of  a 
waiting  time  exceeding  large  t,  which  involves  <,  it  is  desirable  to  appraise  the 
errors  involved.  Two  error  sources  are:  the  systematic  error  (bias)  resulting 
from  fitting  an  incorrect  model,  and  the  effect  of  finite  sample  size  (random 
error).  The  latter  is  the  easiest  to  evaluate,  provided  the  underlying  model  is 
nearly  correct.  The  bias  issue  is  more  difficult  to  deal  with;  one  approach  is 
to  begin  by  fitting  an  elaborate,  multi  parameter  model  and  then  to  check  for  the 
contribution  of  extra  parameters  in  a  prediction  context.  Since  this  paper  takes 
a  non -parametric  approach  to  estimation,  the  bias  resulting  from  an  incorrect 
specification  of  the  service  time  distribution  is  presumably  absent,  although  the 
assumptions  of  stationary  Poisson  arrivals,  and  Independence  are  still  reflected 
in  the  basic  transform  (2.1),  and  hence  in  the  estimates.  Unfortunately,  classi¬ 
cal  procedures  that,  for  example,  avail  themselves  of  the  asymptotic  approxima¬ 
tions  of  properties  of  maximum  likelihood  parameter  estimates  to  estimate  sampling 
errors  are  no  longer  available  in  the  present  environment.  We  have  chosen  in¬ 
stead  to  investigilte  the  performance  of  several  modern  procedures  often  recom¬ 
mended  for  obtaining  standard  errors  of  estimate  and  confidence  intervals:  the 
jackknife  (Quenouille,  Tukey,  Miller,  Hinkley,  and  others)  and  the  bootstrap 
(Efron  and  associates  and  others).  These  methods,  particularly  the  bootstrap, 
which  involves  simulation,  are  computer  intensive,  but  are  the  only  alternatives 
currently  known  for  dealing  with  complex  situations  such  as  that  at  hand. 

The  Jackknife 

■  I  ■  » 

The  above  name  refers  to  a  procedure  originally  introduced  by  Quenouille  for 
bias  reduction  (1956),  and  adapted  by  Tukey  (1958)  to  obtain  approximate  confidence 
intervals.  Suppose  interest  ij^  on  a  parameter  0  (e.g.  our  k,  or  some  function 
thereof)  that  is  estimated  by  0,  using  a  complex  calculation  from  data 
(xi ,X2,x. . ,Xn) ,  just  as  k  is.  The  idea  is  that  of  assessing  variability  by  recom¬ 
puting  K  after  removing  independent  subgroups  of  data  of  equal  size,  and  then 
using  the  recomputed  <  values  to  estimate  a  variance,  which  is  in  turn  applied  to 
state  a  standard  error  or  a  two-sided  confidence  Interval  that  contains  the  true 
0  with  specified  confidence.  A  few  details  follow;  for  more,  see  Efron  (1982) 
and  his  more  recent  work,  or  Mosteller  and  Tukey  (1977).  The  actual  calculation 
involves  splitting  n  into  g  disjoint  groups  of  size  m;  n  =  mg.  Then  calculate 
0fcj),  j  »  1,2,. ..,g:  the  estimate  of  0  that  omits  the  j^  group.  Now  Tukey 
(but  not  Efron)  computes  pseudo-values 

A  A 

yj  »  ge  -  (g-l)0(.j) 
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which  are  then  treated  as  Independent:  use  y  ■  §jk  or  the  point  estimate  of  0, 
and  its  approximate  variance  as 

Sy/g  •  ?  (yj-y)^/(g-i)(g)  •  f 

as  it  turns  out.  Tukey  recommends  referring  y  to  Student's  t  with  g-1  degrees  of 
freedom  to  obtain  confidence  limits.  As  the  name  "jackknife"  Is  Intended  to 
Imply,  the  tool  Is  Inexact  and  a  bit  crude  for  small  samples— just  as  a  true 
jackknife  Is  not  well-adapted  for  delicate  surgery.  It  Is  a  handy  non-pa rametric 
option. 

The  Bootstrap 

In  simplest  form  the  bootstrap  suggests  creating  from  observations  (x^, 

1  *  l,...,n)  the  empirical  distribution 

%<*'  ■  W  ‘<=‘1  '<>  ’ 

Then  one  re-samples  wljh  replacement  from  to  obtain  a  bootstrap  sample  of  size 
n,  from  which  e  (e.g.  k(1))  Is  calculated.  Thls^process  Is  repeated  B  (e.g. 
3Q0-500}  times,  and  the  empirical  distribution,  of  the  resulting  estimates, 
{K(b),  b  *  1,2,.., 8}  Is  employed  as  an  approximation  to  the  sampling  distribution 
of  ic:  the  upper  and  lower  10%  points  of  F"'  approximate  two-sided  80%  confidence 
limits  for  ic,  for  Instance. 


x^  £  X 


x^  >  X 


3.  Transient  Response:  Exponential  Characterizations 

Applications  of  queueing  theory  frequently  require  assessment  of  service  sys¬ 
tem  transient  behavior,  particularly  when  the  system  Is  Initially  Idle  and  heavy 
traffic  conditions  prevail,  for  then  approach  to  steady-state  conditions  may  be 
slow.  Odoni  and  Roth  (1983)  have  reviewed  work  on  the  nature  of  the  approach  to 
steady-state  conditions,  and  have  Investigated  characterizations  of  that  approach 
that  are  simply  stated  and  comprehensible,  being  exponential .  In  this  section 
we  elaborate  upon  the  theme  of  exponential  approach.  Work  by  Abate  and  Whitt 
(1985)  has  similar  alms,  but  proceeds  somewhat  differently. 

3.1  From  Transform  to  Exponential 

Suppose  {X{t),  t£0}  Is  the  state  variable  of  a  process,  and  that  expressions 
for 

(a)  k(s)  -  /  e'*^E[i|/(X(t))|X{0}]dt  »  /  e"^V(t)dt 

*  0  0  * 

(3.1) 

(b)  11m  E[i|;(X(t)|X(0)]  <  « 

t-H» 

are  known.  An  approximate  representation  of 

(c)  4»x(t)  -  EC.|;(X(t))|X(0)] 

Is  desired.  The  proposed  representation  Is  a  simple  exponential  Interpolation  of 
the  form 

ix(t)  -  ij;x(0)e'®^  +  i|,y(«)(l  -e"®^)  ,  (3.2) 

where  8  governs  the  rate  of  approach  to  the  steady-state  value.  The  proposal  Is 

5 


(3.3) 


to  Identify  B  by  minimizing  the  Integrated  squared  error 

00 

A(B)  •  4x(t)  -iPx(0)e'®^+,j,jj(<»)(l  -e"^^)}^w(t)dt  . 

The  weight  function  w(t)  may  be  made  to  focus  upon  values  of  0  appropriate  for 
t-regions  of  interest.  To  optimize  on  B»  study 

00 

»  K  {i|;jj(t)  -i{;jj(0)e"^^  -  -e"®^)}te'^^w(t)dt  (3.4) 

=  0  . 

This  is  equivalent  to 

00  00  00 

/  Ii;„(t)w(t)te"®^dt  =  ij;„(0)  /  e"^®^tw(t)dt  +  iPy(co)  /  (l-e"^^)e'^^tw(t)dt  (3.5) 

0  *  *  0  ^0 

If  w(t)  H  1  the  above  can  be  expressed  in  terms  of  the  derivative  of  ithe  Laplace 
transform  i|)jj(s): 

»  (3-6) 

3 

often  a  transcendental  equation  that  must  be  solved  numerically  for  B.  Clearly, 
if  no  solution  or  multiple  solutions  exist  then  the  simple  uniyersal  form  (3.2) 
is  called  into  question,  but  in  such  cases  weight  functions  are  usefully  invoked. 
If  Interest  centers  on  ultimate  approach  (t  -*■<»)  to  i^y(“)  then  the  smallest  B 
satisfying  (3.6)  is  of  interest.  ^ 

3.2  The  M/G/1  Example. 

As  an  important  Illustration  of  the  above  ideas,  let  X(t)  »  W(t)  be  the 
virtual  waiting  time  in  an  M/G/1  system,  and  consider  ipu(t)  »  E[W(t)lW(0)  =0]. 

Ue  know  that 


^(s)  Poo^^^T 

s 

where 

"oo*"’  •  •sVx'i:i-b(si]  °  s  *  $  ’’ 

b(s)  being  the  transform  of  a  busy  period  duration,  satisfying 

'  b(s)  =  Fjj(s  +  X(1-b(s))  ; 

Fx  is  the  Laplace-Stieltjes  transform  of  the  service  time  or  message  length  dis¬ 
tribution  Fx.  Numerical  solution  of  (3.6)  is  now  possible.  An  approximation  to 
the  (smallest)  value  of  B  governing  the  ultimate  approach  (t  -<•<»)  when  p  +  1,  i.e. 
in  the  heavy-traffic  situation,  is 


3  = 


(3.9) 


In  general,  h'"  (0)  will  Involve  higher  moments  of  the  service  time;  in  case 
X  ■“  Exp(u)  considerable  simplification  occurs  and 


S  »  li(l-p)^  {8(l+p)}'°-® 


(3.10) 


which  will  not,  of  course  agree  exactly  with  formulas  or  numerical  results  given 
previously  by  other  authors,  e.g.  see  Odoni  and  Roth  (1983)  being  a  compromise 
(“Procrustean")  interpolation  of  a  single  exponential  over  the  infinite  t-range. 

In  order  to  focus  upon  a  e-value  relevant  to  ultimate  approach  it  has  been  found 
that  a  weight  procedure  of  the  form  w.^, (t)  »  [l-exp(-6(k)t)]  is  workable  and 
tractable;  start  with  B(0)  *  0  (the  unweighted  procedure  described),  then  intro¬ 
duce  W|^^^  to  determine  and  iterate  to  convergence,  which  occurs  rapidly. 

A  similar  approach  should  apply  to  find  a  B  appropriate  for  small  t. 

The  procedures  described  above  should  be  applicable  in  the  data-driven  context 
of  the  problem  of  Section  1:  a  first  step  is  to  introduce  the  empirical  transform 

A  A  A  A 

and  from  it  proceed  to  b(s),  to  ii<y(s),  and  then  to  6’  Progress  in  this  area, 
and  details  of  the  above  investigations  will  be  reported  elsewhere. 
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